Mean-field study of itinerant ferromagnetism in trapped ultracold Fermi gases: 

Beyond the local density approximation 
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We theoretically investigate the itinerant ferromagnetic transition of a spherically trapped ultra- 
cold Fermi gas with spin imbalance under strongly repulsive interatomic interactions. Our study is 
based on a self-consistent solution of the Hartree-Fock mean-field equations beyond the widely used 
local density approximation. We demonstrate that, while the local density approximation holds in 
the paramagnetic phase, after the ferromagnetic transition it leads to a quantitative discrepancy 
in various thermodynamic quantities even with large atom numbers. We determine the position of 
the phase transition by monitoring the shape change of the free energy curve with increasing the 
polarization at various interaction strengths. 



I. INTRODUCTION 



Recent experimental progress with Feshbach reso- 
nances in ultracold atomic Fermi gases has created op- 
portunities to investigate long-standing many-body prob- 
lems in condensed matter physics. One interesting is- 
sue is the problem of itinerant ferromagnetism in two- 
component (spin-1/2) Fermi gases with repulsive inter- 
actions. The study of itinerant ferromagnetism in con- 
densed matter physics is a fundamental problem which 
has an extensive history, dating back to the basic model 
proposed by Stoner However, the phase transition 
theory of itinerant ferromagnetism is still qualitative. It 
is thought that a Fermi gas with repulsive interactions 
may simulate the Stoner model and therefore could un- 
dergo a ferromagnetic phase transition to a spin-polarized 
state with increased interaction strength [2-6]. This is a 
result of the competition between the repulsive interac- 
tion and the Pauli exclusion principle. The former tends 
to induce polarization, to reduce the interaction energy, 
while the latter prefers a balanced system with equal spin 
populations - and hence a reduced kinetic energy. Above 
a critical interaction strength, the reduced interaction 
energy for a polarized Fermi gas will overcome the gain 
in kinetic energy. Hence, a ferromagnetic phase transi- 
tion should occur. Itinerant ferromagnetism is therefore 
a purely quantum-mechanical effect which occurs when 
the minimum energy is at nonzero magnetization, due to 
the Pauli principle. 

Recently, an_ experimental group reported progress in 
this direction which has attracted intense theoretical 
interest (3, [^4l5|. By using a non-adiabatic field switch 
to the upper branch of a Feshbach resonance with a pos- 
itive s-wave scattering length a s > 0, the experimental- 
ists realized a two-component "repulsive" Fermi gas of 6 Li 
atoms in a harmonic trap, despite the fact that the lower 
branch of the Feshbach resonance is a molecule where the 
fermions attract each other. Initial evidence attributed 
to a ferromagnetic phase transition has been observed, 
suggesting that the transition takes place at k F a s ~ 2.2, 



where k F is the Fermi vector of an ideal gas at the trap 
center. Earlier experiments measuring interaction ener- 
gies [H also provide qualitative evidence in this direction. 

For a trapped Fermi gas, many physical quantities may 
be used to characterize the ferromagnetic phase transi- 
tion. The most direct one is the density profile of each 
component. It was suggested from mean- field theory that 
with a large scattering length, the two Fermi compo- 
nents in a trap should become spatially separated 
This would imply that the majority component stays at 
the trap center and the minority is repelled to the trap 
edge. This spatial inhomogeneity is evidence of a fer- 
romagnetic phase. However, these inhomogeneous den- 
sity profiles or spin domains were not observed in the 
recent experiments, which is away from thermal equilib- 
rium due to the dynamical use of the upper branch of 
Feshbach resonances. For such a non-equilibrium state, 
the onset of phase transition is better characterized by 
a suppression of inelastic three-body collisions, together 
with a minimum in kinetic energy, and a maximum in 
cloud size. From this evidence, a phase transition was 
experimentally identified at a critical scattering length 
of k F a s ~ 2.2 and a temperature of T ~ 0.12Tp. 

On the theoretical side, itinerant ferromagnetism is dif- 
ficult to treat quantitatively. The phase transition occurs 
at large interaction strengths, where fluctuations can be 
huge. Currently, even the critical interaction strength for 
the transition is still under debate. For a homogeneous 
gas, mean-field theory predicts a zero-temperature criti- 
cal value kpa s = tt/2 [l|, while a second-order perturba- 
tive calculation Q predicts a much lower critical coupling 
strength of kp a s ~ 1.054. More accurate quantum Monte 
Carlo simulations recently reported even lower critical 
values k F a s ~ 0.8 [Til. [l4l fl5| ■ For the experimental situ- 
ation with harmonic traps, a local density approximation 
(LDA) is extensively used in order to average over the 
entire trap. At zero temperature, the LDA calculation 
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nite but low temperature causes an increase of the critical 
interaction strength. At the lowest accessible tempera- 
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ture in the recent experiments (T ~ 0.12Tp), we found 
that the LDA predicts k Q F a s ~ 1.93. 

In this paper, rather than addressing the challenging 
problem of calculating the critical interaction strength, 
we instead examine the validity of the widely used lo- 
cal density approximation, by solving the Hartree-Fock 
mean-field equations self-consistently for a trapped im- 
balanced Fermi gas. The LDA is believed to be ap- 
plicable in general cases. However, in extreme condi- 
tions such as a spin-imbalanced superfluid Fermi gas in 
a strongly anisotropic trap, it necessarily breaks down 
due to the strongly distorted superfluid-normal interface, 
which gives a large surface energy [l6l - fl8| . The same sit- 
uation could arise in the spatially phase-separated fer- 
romagnetic phase. In this respect, while the mean-field 
approximation used in this work is not exact, our exam- 
ination of the LDA may provide useful insight for future 
applications of the LDA to more accurate homogeneous 
equations of state of a repulsive Fermi gas. On the other 
hand, beyond-LDA corrections were recently considered 
by LeBlanc and co-workers through the phenomenologi- 
cal inclusion of a surface energy term [6fl . Our study may 
be used to determine accurately the phenomenological 
parameters as an input. 

Our main results may be summarized in the following. 
First, we find a notable disagreement between the self- 
consistent Hartree-Fock solution and the LDA result in 
the strongly interacting ferromagnetic phase where the 
two species are spatially separated, even with large atom 
number up to 10 5 . This discrepancy arises from the non- 
negligible surface energy, which become increasingly im- 
portant at large interaction strength. Secondly, we calcu- 
late various thermodynamic quantities at different tem- 
peratures and spin polarizations. The obtained kinetic 
energy, interaction energy, and the atomic loss rate all 
exhibit a turning point as an indication for the ferromag- 
netic phase transition. We also examine the shape change 
of the free-energy curves as a function of spin polariza- 
tion with increasing the interaction strength, from which 
we are able to determine accurately a mean-field critical 
interaction strength. We find that this does agree with 
the LDA prediction, presumably because surface terms 
play a smaller role at these interaction strengths. 

The rest of the paper is organized as follows. In Sec. 
ILT1 we briefly outline the model and the Hartree-Fock 
mean-field approach. In Sec. IIII1 we present a detailed 
comparison of our results with that obtained from the 
mean-field LDA, as well as a detailed analysis for various 
thermodynamic quantities. Sec. IIVI is devoted to the 
conclusions and further remarks. 



II. METHODS 

We consider an imbalanced Fermi gas with unequal 
spin populations in a spherically trap with an effective 
repulsive interaction, as a minimum model for these ex- 
periments. The system may be described by the Hamil- 



tonian, 

H=J2 j dr¥ a (r)H^ a (r) + U J dA\¥J>^, (1) 

where the pseudo-spin a =t, I denotes different hyperfine 
states and ^ a (r) is the annihilation Fermi field opera- 
tor for the spin-cr state. The single-particle Hamiltonian 
H s a = -fr 2 V 2 /(2m) +V(r) - [i a and V (r) = mw 2 r 2 /2 
is an isotropic harmonic trapping potential with fre- 
quency lj. The effective repulsive interaction strength 
is U = ATrh 2 a s /m > in the lowest Born approxima- 
tion. The total number of fermions and the number 
difference in different hyperfine states are, respectively, 
N = iV-f + N_i and 5N = — N±. We have introduced 
different chemical potentials, /i-^ = /i ± 5[A, to account 
for the number difference or population imbalance. 

We use the Hartree-Fock mean-field (MFA) approx- 
imation to solve the Hamiltonian, which amounts 
to decoupling the interaction term ^I'j^'j^*!'-)- = 

n T (r) ^i>i+n^ (r) (r) (r). Here, n a (r) =< 

a > is the density. The Heisenberg equation of mo- 
tion for ty a (r) then reads 



itldt^a 
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V 2 + T/(r) + [/r^(r)-/i CT 



(2) 



We solve the equations of motion by expanding ^> a (r) = 
Ylj Uj a (r) Cj a exp [— iEja-t], where the field operator Cj a 
annihilates a spin-cr fermion in state j with wave function 
Uj a (r) and energy Ej„. This yields the single-particle 
Schrodinger equation, 
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S7 2 + V(r) + Un w {r)- f i a 



EjaUj a , (3) 



where Uj a (r) are normalized as J d?r\v,j a \ 2 = 1. The 
density distribution n a (r) of spin-cr species can be writ- 
ten as 



(r) 
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(r)i 2 m.), 



(4) 



where / (x) — 1/ [exp (f3x) + 1] is the Fermi distribution 
function at an inverse temperature of j3 = l/(fcsT). The 
chemical potential is determined through the number of 
atoms for the two species: 



J dm a (r) 



= N a . 



(5) 



We self-consistently solve the coupled equations ©- 
([5]). For computational reasons, we introduce a high- 
energy cut-off E c to truncate the summation in the den- 
sity equation For high- lying modes with E > E c , we 
then adopt a semi-classical approximation by assuming 
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that the wavefunction Uj„ is locally a plane wave. We 
refer to refs. [l8[ and (l9| for details. The density profile 
can therefore be written in the form. 



21 + 1, 
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Ki\ 2 f(Kl)+ deU(e),m 



where £ ff = m^fe + [i a — V (r) - Un w (r) / (27r 2 fi 2 ) and 
u ni ( r ) i s * ne ra dial wavefunction. For the given numbers 
of atoms N and SN, temperature T, and s-wave scat- 
tering length a s , our self-consistent iterative procedure 
runs as follows: (i) Start with a Thomas-Fermi profile 
for the density distribution of each species or the pre- 
vious determined density distribution n a (r); (ii) Solve 
the equation ^ for all the radial wave functions u£j(r) 
with energy levels below the cutoff E c ; (iii) Calculate 
the new density profiles, adding both contributions from 
low-lying modes and high-lying modes; (iv) Update the 
chemical potentials according to the number equation ([5]) 
and (v ) finally check the convergence by comparing the 
old and new density profiles. The above steps are re- 
peated until the difference in the old and new density 
profiles decreases to an acceptable level. We note that 
the value of the high-energy cutoff should be carefully 
examined so that the numerical results are independent 
of the choice of E c . 

For the thermodynamic properties, we determine the 
entropy S straightforwardly by 



5=-feE(2l + l)/K)ln/K) ! 



(7) 



where the summation is restricted to the low-lying modes 
below the cutoff E c , as the contribution from the high- 
lying modes is exponentially small. We also calculate the 
total energy E using 



£ E° nl f{E° nl ) + / dr HdeeUie) 

7 ^ TP J J E r 
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-U J dr4nr 2 n t (r) n; (r) + ^ p, a N a 



(8) 



The free energy at finite temperatures is obtained by us- 
ing F = E - TS. 

Experimentally, the kinetic energy E\^ n can be ob- 
tained by measuring the radial width of the ultra-cold 
gas, after free expansion without a trapping potential or 
Feshbach magnetic field. In our numerical work, we cal- 
culate the kinetic energy by subtracting from the total 
energy the potential energy E pot — J drV (r) n (r) and 
the interaction energy E- ia t — U J drn^ (r) (r) , namely, 
Eki n = E — Ep t — -E-int ■ Another observable quantity is 
the three-body loss rate, which, in the weakly interacting 
regime, may be estimated by the expression [20(, 



r = r A' 



j drrif (r) n± (r) n (r) 



(9) 



Here A = kpa s and n (r) = (r) + (r). The prefactor 
Tq may be determined theoretically. While its detailed 
value is not of interest in determining the ferromagnetic 
phase transition, we note that Tq depends on the scatter- 
ing length and Feshbach resonance properties. As this is 
a low-density approximation, finite density effects must 
also be taken into account in interpreting three-body loss 
data as a guide to atomic density and polarization. 



III. RESULTS AND DISCUSSIONS 

In the numerical calculations, we use the "trap units" 
h = m = uj = k^ = \. The length, energy and temper- 

1 /2 

ature are then measured in units of = [h/ (raw)] , 
Hlo and Hco/kB, respectively. Furthermore, we will scale 
our results using characteristic quantities obtained from a 
zero-temperature ideal balanced Fermi gas in an isotropic 
harmonic trap. These are the Fermi energy Ep = 

1/3 

(3iV) Huj, Fermi temperature Tp = Ep/ks, Thomas- 
Fermi radius rpF = (24iV) 1 ^ 6 a^o and the peak density 
npF = (6iV) 1//2 /3ir 2 a^. In most cases, we set the total 
atom number = I0 5 and temperature T — 0.12Tp, 
comparable to the parameters used in the recent MIT 
experiment. To remove complications due to spin degen- 
eracy, we also set a small polarization p = (Nf — N±)/N = 
0.10. We use a cut-off energy E c = 85Huj and solve the 
Hartree-Fock mean-field equation within the subspace 
"■max = 80 and / ma x = 120, which contains all the en- 
ergy levels with energy E < E c . 

In the following, we first compare the Hartree-Fock 
mean-field density profiles with the corresponding LDA 
predictions to examine the validity of LDA. We then in- 
vestigate in detail some thermodynamic properties and 
determine the critical interaction strength for emergence 
of the ferromagnetic state. 



A. Hartree-Fock MF versus LDA 

In Fig. [TJ we present the density profile of each spin 
component and the total density profile of the imbalanced 
Fermi gas, calculated from Hartree-Fock MF (solid lines) 
and LDA (dashed lines) methods with total atom num- 
ber N = 10 5 and polarization p = (iV t - N±)/N = 0.10 
at T = Q.12Tp. When the interaction strength is rela- 
tively weak (e.g., before the ferromagnetic phase transi- 
tion, k F a s — 2.0, left column), the LDA prediction agrees 
fairly well with that from the Hartree-Fock MF theory. 
Similar excellent agreement is also found in all the phys- 
ical quantities that we have considered, including the 
chemical potential, energy, and the atom loss rate (see, 
for example, Fig. [3]). However, when the repulsive inter- 
action is stronger than a critical value (e.g., k F a s = 4.0, 
right column), there is a significant discrepancy. In this 
parameter space, we observe a spatial phase separation 
in the density profiles, with majority spin-up atoms oc- 
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Figure 1: (Color online) Spin-up, spin-down and total density 
profiles of an imbalanced Fermi gas at k° F a a — 2.0 (left col- 
umn) and kpd s = 4.0 (right column). The solid and dashed 
lines refer to the Hartree-Fock MF and LDA results, respec- 
tively. The inset in the upper, right-column plot highlights 
the spin-up density profile at the trap edge. The spin polar- 
ization is p = 0.10 and the number of atoms is TV = 10 5 . 



cupying the inner core and minority spin-down atoms 
repelled to the edge of the trap. This density separation 
is revealed clearly by both the LDA calculations and the 
Hartree-Fock MF calculations. However, LDA predicts 
a sudden change of the density profile of each spin com- 
ponent at r ~ 0.7rTF, due to the neglect of the spatial 
derivatives of the densities. At this point, the LDA neces- 
sarily fails and one has to resort to the more reliable self- 
consistent Hartree-Fock theory. Another notable discrep- 
ancy between LDA and Hartree-Fock MF comes from the 
trap edge r ~ ttf, where the Hartree-Fock MF theory 
gives a much larger spin-up density than the LDA (see, 
for example, the inset in the right, upper plot). This is a 
well-known breakdown of the LDA method due to small 
densities occurring at the trap edge. Because of this fail- 
ure, as we shall see, the LDA predicts a much smaller 
atom loss rate than the Hartree-Fock method. We note 
also that, for r ~ ttf the density of both two species be- 
comes very low. There is no spatial phase separation, as 
a result of much weaker effective repulsive interactions. 

We have also checked a smaller Fermi cloud with 
N = 10 3 at the same temperature and polarization. As 
anticipated, the difference between Hartree-Fock MF and 
LDA becomes even more significant. 
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Figure 2: (Color online) Magnetization profile or local spin 
polarization at different interaction strengths: k^a s = 1.8 
(dash-dotted line), 2.7 (dashed line), and 4.0 (solid line). 
These curves are calculated using the Hartree-Fock MF the- 
ory for a total number of atoms N = 10 and a polarization 
p = 0.10 at T = 0.12T F . 



B. Magnetization profile 

As shown above, to reliably characterize the spin den- 
sity or spin domains in the ferromagnetic phase, one must 
include corrections beyond-LDA by considering surface 
energy. This was done in the ref. [6| by phenomenologi- 
cally including a magnetization gradient term. 

In Fig. 2, we report the Hartree-Fock MF results for 
the magnetization profile, m (r) = (n-j- (r) —n^(r))/n (r). 
Before the transition (k F a s — 1.8), the magnetization 
is nonzero but small. When the repulsive interaction is 
larger than the critical value (k F a s = 2.7 and 4.0), we 
find a fully magnetized core, whose size increases with 
increasing the interaction strength. As the strength in- 
creases, another nearly fully magnetized domain forms 
in the vicinity of the trap edge with m (r) w — 1 . These 
observations are in qualitative agreement with the theo- 
retical predictions reported in ref. [6j. 



C. Thermodynamic properties 

We now turn to thermodynamic properties of the im- 
balanced Fermi gas with repulsive interactions, and com- 
pare our results with the experimental measurements 
where available. In Fig. ([3]), we show the atom loss rate 
(a), interaction energy (b), potential energy (c), kinetic 
energy (d) and chemical potentials (e and f ) . It is readily 
seen that LDA and Hartree-Fock MF calculations agree 
extremely well with each other for relatively weak interac- 
tions below the critical interaction strength (k F a s ~ 2.0), 
as mentioned earlier. However, the agreement tends to 
be worse as the interaction strength becomes stronger. 

The most significant discrepancy occurs in the atom 
loss rate of three-body inelastic collisions. In the Hartree- 
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Figure 3: (Color online) Atom loss rate (a), interaction en- 
ergy (b), potential energy (c), kinetic energy (d), and chemi- 
cal potentials (e and f) as a function of interaction strength. 
The Hartree-Fock MF results (solid lines) are compared with 
the LDA predictions (dashed lines). The vertical dashed line 
marks the ferromagnetic phase transition. Here, N = 10 , 
p = 0.10, and T = 0.127>. 



Fock MF calculation, the loss rate increases monoton- 
ically with increasing interaction strength. While in 
the LDA calculations, it increases first but then de- 
creases when the interaction strength is above the crit- 
ical value, as illustrated in Fig. ([3^). This discrep- 
ancy lies in the different prediction of the density tail 
from the two methods. It is clear that the loss rate 
is proportional to the integrated density profile overlap 
rif (r) nj, (r)(n-f (r) + (r)) over the entire space. For 
strong repulsion, the spin-up and spin-down density pro- 
files are essentially not overlapping (i.e., spatially phase 
separated), apart from the small regions at the phase- 
separation boundary, where the Hartree-Fock MF calcu- 
lation shows a residue overlap, and at the edge of the 
trap, where due to the small density, the effective repul- 
sion cannot sustain the phase separation further. The 
atom loss therefore arises mostly at the phase-separation 
boundary and at the trap edge. In these two regions, the 
LDA fails to give a quantitative density distribution and 
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Figure 4: (Color online) Atom loss rate (a), interaction en- 
ergy (b), potential energy (c), and kinetic energy (d) at dif- 
ferent temperatures T = 0.22T F (solid lines) and T = 0.557> 
(dashed line) , with a total number of atoms N = 10 4 and an 
imbalance ratio p = 0.10. 



hence a reliable prediction for the atom loss rate. 

However, it should be noted that, the qualitative ex- 
perimental observation [2j, the suppression of the atom 
loss rate at large repulsion, agrees with the LDA pre- 
dictions, but is in contradiction with our more accurate 
Hartree-Fock MF calculations. We attribute the discrep- 
ancy to the applicability of equation ([5]), in which the 
strong dependence of the loss rate on the interaction 
strength (i.e., scaling as (kpa s ) 6 ) significantly overesti- 
mates the loss rate above the phase transition. This 
dependence, however, is only valid in the weakly inter- 
acting regime (kpa s <C 1). A non-perturbative theo- 
retical expression for the loss rate is therefore needed for 
large interaction strength, with which we expect that the 
Hartree-Fock MF theory would then predict a suppres- 
sion of the atom loss rate. 

There are very clear quantitative differences in the in- 
teraction energy, potential energy and kinetic energy in 
the ferromagnetic phase predicted by the two methods 
(see Figs. (|3Jd), ((3t), and (J3J1)). These differences are 
largely significant once the ferromagnetic domains are 
formed, and presumably are due to errors in the LDA 
treatment of the domain boundaries. Qualitatively, all 
the results from the two methods indicate a ferromag- 
netic phase transition at about k° F a s w 2.0. The inter- 
action energy and potential energy show a maximum at 
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the transition, while the kinetic energy exhibits a min- 
imum. Quantitatively, our mean-field prediction of the 
critical interaction strength is (k F a s ) c — 1.93. This is 
in good agreement with the MIT experimental finding of 
(k F a s ) c ~ 2.2 at the same temperature T — 0.12TV. 

However, this agreement may simply be a coincidence. 
Monte-Carlo methods that take fluctuations into account 
predict a lower critical strength. Mean-field theory is not 
quantitatively reliable at these large interaction strengths 
with k F a s > 1. A possible explanation for the appar- 
ent coincidence is that the experimentally determined 
critical strength may not accurately correspond to the 
true, equilibrium critical strength. The recent measure- 
ments were carried out dynamically. In this type of non- 
equilibrium state, critical slowing-down and hysteresis- 
like effects will prevent ferromagnetic domains from im- 
mediately forming, and may push the effective critical 
interaction strength to a higher value. 

The critical interaction strength also depends crucially 
on temperature. In Fig. ([4}, we graph the temperature 
dependence of the atom loss rate and energies. At higher 
temperatures, T — 0.22XV, the mean-field critical inter- 
action strength increases to (k F a s ) ~ 2.4,. This is much 
smaller than the experimental result of (k F a s ) c ~ 4.2, at 
the same temperature. For an even higher temperature, 
T _± 0.557V, there is no non-monotonic behavior in our 
mean-field results of the atomic loss and energies and, 
therefore no indication for the phase transition anymore. 
This is in agreement with experimental observation. 

We finally consider the chemical potentials as shown 
in Figs. ((3^) and ((3f). The chemical potential increases 
with interaction strength and nearly saturates in the fer- 
romagnetic phase. The same saturation was qualitatively 
observed in the MIT experiment jjjj. 



D. Critical interaction strength from the shape of 
free energy curve 
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Figure 5: (Color online) (a) Free energy as a function of spin 
polarization p for different interaction strengths: k F a B = 1.94, 
2.00, 2.10, 2.20, 2.40 and 3.0. (b) Free energy vs interaction 
strength k^a s for different spin imbalances: p — 0.05, 0.10 
and 0.20. Here, N = 10 4 and T = 0.122V. 



At a finite temperature, we can accurately determine 
the critical interaction strength from the free energy, 
which should exhibit a minimum at non-zero magnetiza- 
tion. We present in Fig. ([5^) the free energy as a function 
of spin polarization at different interaction strengths. For 
a weak interaction strength, for example k F a — 1.92, the 
free energy increases monotonically with spin polariza- 
tion. Since the free energy should be symmetric (same) 
for ±p, we simply refer this monotonic dispersion to as 
"U" shape. When the interaction is large enough, the 
state with a nonzero spin polarization becomes energeti- 
cally favored (see, e.g., k F a s = 2.4). Thus, the dispersion 
is changed to a "W" shape. As shown in Fig. (^S^,), the 
position of the minimum shifts to large spin polarization 
with increasing interaction strength. 

To determine the interaction strength accurately for 
this shape change, we calculate the free energy F p at 
three small spin polarizations: p = 0, 0.02 and 0.05 near 
k F a s = 2.0. In Fig. (03), and show the free energy with 



respect to the non-polarization value F p= $. We identify 
a crossover point at k F a s — 1.93, above (below) which 
the state with spin polarization is lower (higher) in free 
energy. This is the critical interaction strength. For 
these atom numbers, our mean-field critical strength is 
the same as the LDA prediction for a trapped repulsive 
Fermi gas at the same temperature. 



IV. CONCLUSIONS 

In conclusion, we have performed a self-consistent 
Hartree-Fock mean-field study of itinerant ferromag- 
netism for a harmonically trapped Fermi gas with 
strongly repulsive interactions. Our results for the den- 
sity profiles, equations of state, and atom loss rate have 
been used to examine the validity of the widely used lo- 
cal density approximation. We have found that the local 
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density approximation works fairly well below the ferro- 
magnetic phase transition, and gives the same prediction 
for the critical coupling strength. However, in the fer- 
romagnetic or spatially inhomogeneous phase, it gives 
quantitatively different predictions for both the equa- 
tions of state and the atom loss rate, even for numbers of 
atoms as large as N = 10 5 . This discrepancy between 
the local density approximation and the more precise 
Hartree-Fock mean-field predictions is apparently due to 
the breakdown of the local density approximation at both 
the phase-separation boundary and at the edge of the 
trap. 

This failure of the local density approximation is likely 
to remain an important issue even if we apply the lo- 
cal density approximation to the more accurate quantum 
Monte Carlo data obtained from a homogeneous ferro- 
magnetic Fermi system. For this reason, it would be very 
useful to perform an inhomogeneous quantum simulation 
in the presence of a harmonic trap, in order to extend the 
present work beyond the mean-field approximation. 

Our study may be useful for calculating the stiffness 
in the magnetization gradient term, which is also beyond 
the LDA. The stiffness has been computed phenomeno- 



logically in ref. Q using a Landau-Ginzburg expansion 
for the excess energy. For the normal-superfluid interface 
of a population unbalanced Fermi gas in the unitarity 
limit, it was shown that this expansion of gradient terms 
leads to about a factor-of-two error at zero temperature, 
compared with the full mean-field Bogoliubov-de Gennes 
calculation [2lJ]. We thus anticipate that our Hartree- 
Fock mean-field calculation will determine a similar stiff- 
ness as derived in [(|, with a similar order of magnitude. 
However, a complete determination of the stiffness term 
is numerically challenging [2l|, especially at finite tem- 
peratures. 
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